{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/matplotlib/__init__.py:846: MatplotlibDeprecationWarning: \n",
      "The text.latex.unicode rcparam was deprecated in Matplotlib 2.2 and will be removed in 3.1.\n",
      "  \"2.2\", name=key, obj_type=\"rcparam\", addendum=addendum)\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/matplotlib/__init__.py:855: MatplotlibDeprecationWarning: \n",
      "examples.directory is deprecated; in the future, examples will be found relative to the 'datapath' directory.\n",
      "  \"found relative to the 'datapath' directory.\".format(key))\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/matplotlib/__init__.py:947: UserWarning: could not find rc file; returning defaults\n",
      "  warnings.warn(message)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as npl\n",
    "import scipy.linalg as spl\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12.1 Least Squares Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.816496580927726"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "matrix([[-0.33333333, -0.66666667,  0.33333333]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A,b = np.matrix([[2,0],[-1,1],[0,2]]), np.array([1,0,-1])\n",
    "xhat = np.array([1/3,-1/3])\n",
    "rhat = np.matmul(A,xhat) - b\n",
    "npl.norm(rhat)\n",
    "x = np.array([1/2,-1/2])\n",
    "r = np.matmul(A,x) - b\n",
    "npl.norm(r)\n",
    "rhat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12.2 Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.33333333, -0.33333333]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.33333333, -0.33333333]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "matrix([[-2.22044605e-16,  2.22044605e-16]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.matmul(npl.inv(A.transpose()*A)*A.transpose(),b)  \n",
    "np.matmul(npl.pinv(A),b)\n",
    "np.matmul(np.matmul(A.transpose(),A),xhat) - np.matmul(A.transpose(),b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[2.22044605e-16]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z = np.array([-1.1,2.3])\n",
    "# np.shape(np.matmul(A,z))\n",
    "# np.shape(rhat)\n",
    "np.matmul(A,z) @ rhat.transpose() #J: (A*z)'*rhat; auto columning in Julia means we have to reverse the transpose order "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12.3 Solving Least Squares Problems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "A,b = np.random.randn(100,20),np.random.randn(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:1: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "6.185158925429641e-16"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "4.799638480896434e-16"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "5.453688670685701e-16"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x1 = npl.lstsq(A,b)[0] #lstsq gives the solution, residuals and rank, need to choose [0], \n",
    "                       #in Julia the backslash operator auto uses lstsq based on if given matrix is overdetermined \n",
    "x2 = np.matmul(npl.inv(np.matmul(A.transpose(),A)),np.matmul(A.transpose(),b))\n",
    "x3 = np.matmul(npl.pinv(A),b)\n",
    "Q,R = npl.qr(A)\n",
    "x4 = npl.solve(R,(np.matmul(Q.transpose(),b)))\n",
    "\n",
    "npl.norm(x1-x2)  \n",
    "npl.norm(x2-x3)\n",
    "npl.norm(x3-x4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Complexity: \n",
    "1. doubling m doubles compute time\n",
    "2. doubling n quadruples compute time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:1: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "81.1 ms ± 8 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "m,n = 2000, 500\n",
    "A,b = np.random.randn(m,n), np.random.randn(m)\n",
    "%timeit npl.lstsq(A,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:1: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "124 ms ± 14.8 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "m,n = 4000, 500\n",
    "A,b = np.random.randn(m,n), np.random.randn(m)\n",
    "%timeit npl.lstsq(A,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:1: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "311 ms ± 43.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "m,n = 2000, 1000\n",
    "A,b = np.random.randn(m,n), np.random.randn(m)\n",
    "%timeit npl.lstsq(A,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:2: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:3: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3.614235711934445e-16"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A,B = np.random.randn(1000,100), np.random.randn(1000,10)\n",
    "X = npl.lstsq(A,B)[0]\n",
    "x3 = npl.lstsq(A,B[:,3])[0]\n",
    "npl.norm(X[:,3]-x3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12.4 Examples "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:4: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  after removing the cwd from sys.path.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([  62.07662454,   99.98500403, 1442.83746254])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#advertising budget\n",
    "R = np.matrix([[.97,1.86,.41],[1.23,2.18,.53],[.8,1.24,.62],[1.29,.98,.51],[1.1,1.23,.69],[.67,.34,.54],[.87,.26,.62],[1.1,.16,.48],[1.92,.22,.71],[1.29,.12,.62]])\n",
    "m,n = np.shape(R)\n",
    "vdes = 1e3 * np.ones(m)\n",
    "s = npl.lstsq(R,vdes)[0]\n",
    "s\n",
    "#will be continued in 16 re: how to add a constraint like a total budget "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "132.63819026326522"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rms = lambda x: np.sqrt(np.mean(np.square(x))) #not a built in numpy function\n",
    "rms(np.matmul(R,s) - vdes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 625)"
      ]
     },
     "execution_count": 143,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#10 lamps each with an x,y position and a height above the floor\n",
    "#goal: adjust light intensity for each lamp to get uniform lighting at the floor level\n",
    "lamps = np.matrix([[4.1,20.4,4],[14.1,21.3,3.5],[22.6,17.1,6],[5.5,12.3,4.0],[12.2,9.7,4.0],[15.3,13.8,6],[21.3,10.5,5.5],[3.9,3.3,5.0],[13.1,4.3,5.0],[20.3,4.2,4.5]])\n",
    "n = np.shape(lamps)[0]\n",
    "N = int(max(np.ravel(lamps))+2.4) #grid size\n",
    "m = N*N #pixel density\n",
    "n,m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [],
   "source": [
    "#an m x 2 matrix with coordinates of pixel centers \n",
    "pixels = np.hstack((np.vstack([np.vstack(np.arange(.5,N,1)) for i in range(N)]),np.vstack([np.vstack(np.full(N,.5+i)) for i in range(N)])))\n",
    "A = np.zeros((m,n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:5: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \"\"\"\n"
     ]
    }
   ],
   "source": [
    "for i in range(m):\n",
    "    for j in range(n):\n",
    "        A[i,j] = 1.0/np.square(npl.norm(np.hstack([pixels[i,:],[0]]) - lamps[j,:]))\n",
    "A = m/sum(A) * A\n",
    "x = npl.lstsq(A,np.ones((m,1)))[0]\n",
    "rmsLS = rms(np.matmul(A,x) - 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "rms(np.matmul(A,np.ones((n,1))) - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  7.,  13.,  35.,  58.,  78., 168., 156.,  69.,  35.,   6.]),\n",
       " array([0.53478587, 0.61503502, 0.69528417, 0.77553332, 0.85578247,\n",
       "        0.93603162, 1.01628077, 1.09652992, 1.17677907, 1.25702821,\n",
       "        1.33727736]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 157,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEZdJREFUeJzt3X+w5XVdx/HnKzY0LQPdi+EutFiLujo40pWhNEPpB4jDYmqzjD/WonY0MvvhBMRMODVMazVaTubMhhtrYyiRCoVWRBiTCnqR379kQ4Irq3sNwcoGRd/9cb5btzt399x7fuw5fHo+Zu7c8/18P+d8X3Pv3dd+7/d+v9+TqkKS1K7vmHQASdJ4WfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxq2ZdACAtWvX1oYNGyYdQ5IeV2644YavVNVMv3lTUfQbNmxgbm5u0jEk6XElyb+uZJ6HbiSpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXFTcWWsNM02nHvlRLZ73/bTJrJdtcc9eklqnEUvSY2z6CWpcRa9JDXOopekxvUt+iQ7k+xNctuS8bckuTvJ7Ul+b9H4eUl2d+t+ahyhJUkrt5LTKy8G/hh4/76BJC8FNgPHVdWjSY7oxjcBW4DnAs8A/iHJsVX1rVEHlyStTN89+qq6FnhoyfCbge1V9Wg3Z283vhn4YFU9WlVfAHYDJ4wwryRplQY9Rn8s8KNJrk/yT0le2I2vAx5YNG++G5MkTcigV8auAQ4HTgReCFya5JlAlplby71Akm3ANoCjjz56wBiSpH4G3aOfBz5cPZ8Bvg2s7caPWjRvPfDgci9QVTuqaraqZmdm+r6JuSRpQIMW/UeBlwEkORY4FPgKcAWwJckTkhwDbAQ+M4qgkqTB9D10k+QS4CRgbZJ54AJgJ7CzO+XyG8DWqirg9iSXAncAjwFne8aNJE1W36KvqjP3s+p1+5l/IXDhMKEkSaPjbYqlKeXtkTUq3gJBkhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxvUt+iQ7k+zt3k1q6bq3Jakka7vlJHl3kt1Jbkly/DhCS5JWbiV79BcDpywdTHIU8BPA/YuGT6X3PrEbgW3Ae4ePKEkaRt+ir6prgYeWWfUu4DeAWjS2GXh/9VwHHJbkyJEklSQNZKBj9ElOB75YVTcvWbUOeGDR8nw3ttxrbEsyl2RuYWFhkBiSpBVYddEneRJwPvBby61eZqyWGaOqdlTVbFXNzszMrDaGJGmFBnlz8B8AjgFuTgKwHvhckhPo7cEftWjueuDBYUNKkga36j36qrq1qo6oqg1VtYFeuR9fVV8CrgDe0J19cyLwSFXtGW1kSdJqrOT0ykuATwPPSjKf5KwDTP8YcC+wG/hT4BdHklKSNLC+h26q6sw+6zcselzA2cPHkiSNilfGSlLjLHpJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1LiVvPHIziR7k9y2aOz3k9yV5JYkH0ly2KJ15yXZneTuJD81ruCSpJVZyR79xcApS8auAp5XVccBnwfOA0iyCdgCPLd7zp8kOWRkaSVJq9a36KvqWuChJWN/X1WPdYvX0XsTcIDNwAer6tGq+gK9txQ8YYR5JUmrNIpj9D8HfLx7vA54YNG6+W5MkjQhQxV9kvOBx4AP7BtaZlrt57nbkswlmVtYWBgmhiTpAAYu+iRbgVcAr+3eFBx6e/BHLZq2HnhwuedX1Y6qmq2q2ZmZmUFjSJL6GKjok5wCnAOcXlVfX7TqCmBLkickOQbYCHxm+JiSpEGt6TchySXAScDaJPPABfTOsnkCcFUSgOuq6k1VdXuSS4E76B3SObuqvjWu8JKk/voWfVWduczw+w4w/0LgwmFCSZJGxytjJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mN61v0SXYm2ZvktkVjT01yVZJ7us+Hd+NJ8u4ku5PckuT4cYaXJPW3kj36i4FTloydC1xdVRuBq7tlgFPpvU/sRmAb8N7RxJQkDapv0VfVtcBDS4Y3A7u6x7uAMxaNv796rgMOS3LkqMJKklZv0GP0T6+qPQDd5yO68XXAA4vmzXdjkqQJGfUfY7PMWC07MdmWZC7J3MLCwohjSJL2GbTov7zvkEz3eW83Pg8ctWjeeuDB5V6gqnZU1WxVzc7MzAwYQ5LUz6BFfwWwtXu8Fbh80fgburNvTgQe2XeIR5I0GWv6TUhyCXASsDbJPHABsB24NMlZwP3Aa7rpHwNeDuwGvg787BgyS5JWoW/RV9WZ+1l18jJzCzh72FCSpNHxylhJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuOGKvokv5rk9iS3JbkkyROTHJPk+iT3JPlQkkNHFVaStHoDF32SdcAvA7NV9TzgEGAL8A7gXVW1EfgqcNYogkqSBjPsoZs1wHclWQM8CdgDvAy4rFu/CzhjyG1IkoYwcNFX1ReBP6D35uB7gEeAG4CHq+qxbto8sG655yfZlmQuydzCwsKgMSRJfQxz6OZwYDNwDPAM4MnAqctMreWeX1U7qmq2qmZnZmYGjSFJ6mOYQzc/Dnyhqhaq6pvAh4EfAQ7rDuUArAceHDKjJGkIa/pP2a/7gROTPAn4L+BkYA64Bng18EFgK3D5sCGlDedeOekI0uPWMMfor6f3R9fPAbd2r7UDOAf4tSS7gacB7xtBTknSgIbZo6eqLgAuWDJ8L3DCMK8rSRodr4yVpMZZ9JLUOItekhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNW6oC6YktWeSt5u4b/tpE9t2y9yjl6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY0bquiTHJbksiR3JbkzyQ8neWqSq5Lc030+fFRhJUmrN+we/R8Bf1tVzwaeD9wJnAtcXVUbgau7ZUnShAxc9EmeAryE7q0Cq+obVfUwsBnY1U3bBZwxbEhJ0uCG2aN/JrAA/FmSG5NclOTJwNOrag9A9/mIEeSUJA1omKJfAxwPvLeqXgD8J6s4TJNkW5K5JHMLCwtDxJAkHcgwRT8PzFfV9d3yZfSK/8tJjgToPu9d7slVtaOqZqtqdmZmZogYkqQDGfimZlX1pSQPJHlWVd0NnAzc0X1sBbZ3ny8fSVJNhUne8ErSYIa9e+VbgA8kORS4F/hZer8lXJrkLOB+4DVDbkOSNIShir6qbgJml1l18jCvK0kaHa+MlaTGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1buiiT3JIkhuT/E23fEyS65Pck+RD3btPSZImZBR79G8F7ly0/A7gXVW1EfgqcNYItiFJGtBQRZ9kPXAacFG3HOBlwGXdlF3AGcNsQ5I0nGH36P8Q+A3g293y04CHq+qxbnkeWLfcE5NsSzKXZG5hYWHIGJKk/Rm46JO8AthbVTcsHl5mai33/KraUVWzVTU7MzMzaAxJUh9rhnjui4DTk7wceCLwFHp7+IclWdPt1a8HHhw+piRpUAPv0VfVeVW1vqo2AFuAf6yq1wLXAK/upm0FLh86pSRpYOM4j/4c4NeS7KZ3zP59Y9iGJGmFhjl08z+q6hPAJ7rH9wInjOJ1JUnD88pYSWrcSPbodXBtOPfKSUeQxmJSP9v3bT9tIts9WNyjl6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJatww7xl7VJJrktyZ5PYkb+3Gn5rkqiT3dJ8PH11cSdJqDbNH/xjw61X1HOBE4Owkm4BzgauraiNwdbcsSZqQYd4zdk9Vfa57/O/AncA6YDOwq5u2Czhj2JCSpMGN5Bh9kg3AC4DrgadX1R7o/WcAHDGKbUiSBjN00Sf5buCvgF+pqq+t4nnbkswlmVtYWBg2hiRpP4Yq+iTfSa/kP1BVH+6Gv5zkyG79kcDe5Z5bVTuqaraqZmdmZoaJIUk6gGHOugnwPuDOqnrnolVXAFu7x1uBywePJ0ka1jBvDv4i4PXArUlu6sZ+E9gOXJrkLOB+4DXDRZQkDWPgoq+qfwayn9UnD/q6kqTR8spYSWqcRS9JjbPoJalxFr0kNW6Ys27+39tw7pWTjiBJfblHL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOItekhr3uD+90lMcJenAHvdFL0nDmuQO433bTxv7Njx0I0mNs+glqXEWvSQ1bmxFn+SUJHcn2Z3k3HFtR5J0YGMp+iSHAO8BTgU2AWcm2TSObUmSDmxce/QnALur6t6q+gbwQWDzmLYlSTqAcRX9OuCBRcvz3Zgk6SAb13n0y71peP2fCck2YFu3+B9J7gbWAl8ZU6ZhmGt1zLU65lqdpnLlHUNt8/tXMmlcRT8PHLVoeT3w4OIJVbUD2LF4LMlcVc2OKdPAzLU65lodc62OuVZvXIduPgtsTHJMkkOBLcAVY9qWJOkAxrJHX1WPJfkl4O+AQ4CdVXX7OLYlSTqwsd3rpqo+BnxslU/b0X/KRJhrdcy1OuZaHXOtUqqq/yxJ0uOWt0CQpMZNpOj73R4hyRuTLCS5qfv4+WnI1c35mSR3JLk9yV9MQ64k71r0tfp8koenJNfRSa5JcmOSW5K8fEpyfX+Sq7tMn0iy/iBk2plkb5Lb9rM+Sd7dZb4lyfHjzrTCXM9O8ukkjyZ528HItMJcr+2+Trck+VSS509Jrs1dppuSzCV58cHI1VdVHdQPen+c/RfgmcChwM3ApiVz3gj88RTm2gjcCBzeLR8xDbmWzH8LvT9+TzwXvWOWb+4ebwLum5Jcfwls7R6/DPjzg5DrJcDxwG37Wf9y4OP0rkE5Ebh+3JlWmOsI4IXAhcDbDkamFeb6kUX/Dk+doq/Xd/O/h8SPA+46WF+zA31MYo9+Wm+PsJJcvwC8p6q+ClBVe6ck12JnApdMSa4CntI9/l6WXEsxwVybgKu7x9css37kqupa4KEDTNkMvL96rgMOS3LkpHNV1d6q+izwzXFnWbLdfrk+te/fIXAdvWt1piHXf1TX8sCTWXKh6KRMouhXenuEV3W/Al2W5Khl1k8i17HAsUk+meS6JKdMSS6gd0gCOAb4xynJ9XbgdUnm6Z2B9ZYpyXUz8Kru8SuB70nytIOQ7UC8bcjgzqL329BUSPLKJHcBVwI/N+k8MJmi73t7BOCvgQ1VdRzwD8CusadaWa419A7fnERvz/miJIdNQa59tgCXVdW3xphnn5XkOhO4uKrW0zs08edJxv0zt5JcbwN+LMmNwI8BXwQeG3OuflbzfVYnyUvpFf05k86yT1V9pKqeDZwB/M6k88Bkin4lt0f4t6p6tFv8U+CHpiFXN+fyqvpmVX0BuJte8U861z5bODiHbWBluc4CLgWoqk8DT6R3P5CJ5qqqB6vqp6vqBcD53dgjY87Vz2q+zwKSHAdcBGyuqn+bdJ6lusM8P5Bk3D/zfU2i6PveHmHJscnTgTunIRfwUeClXca19A7l3DsFuUjyLOBw4NNjzrOaXPcDJ3f5nkOv6BcmnSvJ2kW/WZwH7BxzppW4AnhDd/bNicAjVbVn0qGmVZKjgQ8Dr6+qz086zz5JfjBJusfH0zshYPL/CU3iL8D0fo3/PL2zI87vxn4bOL17/LvA7fSOpV4DPHtKcgV4J3AHcCuwZRpydctvB7ZP2fdxE/DJ7vt4E/CTU5Lr1cA93ZyLgCcchEyXAHvo/VFznt5vO28C3rToZ+s9XeZbgdmD9LXql+v7uvGvAQ93j58yBbkuAr7a/VzdBMxNydfrnK67bqK30/Xig5Gr34dXxkpS47wyVpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktS4/wYdxa7EJjuuDgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(np.matmul(A,x)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 10.,  25.,  51.,  80.,  87., 113.,  94.,  99.,  50.,  16.]),\n",
       " array([ 4.35112509,  5.38933026,  6.42753542,  7.46574059,  8.50394576,\n",
       "         9.54215092, 10.58035609, 11.61856126, 12.65676642, 13.69497159,\n",
       "        14.73317676]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 158,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADXZJREFUeJzt3X+o3fV9x/Hna2ZidRO1ubo00V07gp3IOuUitkIZpgVtxOSPCpZuC10g/3StazdqXGH+NyIrsx0bjqDWjImrZB2R6joltZTBKkvU+ivtEqyLV1NzS2u3tQMb+t4f99txiTcmOd9z8r33k+cD5Jzv93zP+b6/Gp/3m++959xUFZKkdv3S0ANIkibL0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDVuxdADAKxcubKmp6eHHkOSlpW9e/f+oKqmjrfdkgj99PQ0e/bsGXoMSVpWkvzniWznpRtJapyhl6TGGXpJapyhl6TGGXpJapyhl6TGGXpJapyhl6TGGXpJatySeGestFRNb31ksH2/vG39YPtWWzyjl6TGGXpJapyhl6TGGXpJapyhl6TGGXpJapyhl6TGGXpJapyhl6TGGXpJapyhl6TGGXpJapyhl6TGGXpJapyhl6TGGXpJapyhl6TGGXpJapyhl6TGHTf0Se5LcjjJ8wvWXZDk8ST7u9vzu/VJ8ldJDiR5NslVkxxeknR8J3JGfz9w/VHrtgK7q2otsLtbBrgBWNv9swW4ezxjSpJGddzQV9U3gR8etXoDsKO7vwPYuGD939W8bwHnJVk1rmElSSdv1Gv0F1XVIYDu9sJu/WrglQXbzXbrJEkDGfc3Y7PIulp0w2RLkj1J9szNzY15DEnSL6wY8XmvJ1lVVYe6SzOHu/WzwMULtlsDvLbYC1TVdmA7wMzMzKJfDCSdetNbHxls3y9vWz/Yvls26hn9w8Cm7v4mYNeC9b/f/fTNNcCPf3GJR5I0jOOe0Sd5EPgdYGWSWeAOYBvwUJLNwEHg5m7zR4EPAweAnwIfn8DMkqSTcNzQV9VHj/HQukW2LeATfYeSJI2P74yVpMYZeklqnKGXpMYZeklqnKGXpMYZeklq3KjvjJU0YUO+Q1Vt8Yxekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhrnZ91oWfBzX6TReUYvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY3rFfokn07yQpLnkzyY5KwklyZ5Msn+JF9Ocua4hpUknbyRQ59kNfApYKaqrgDOAG4B7gTuqqq1wI+AzeMYVJI0mr6XblYA70iyAjgbOARcB+zsHt8BbOy5D0lSDyN/qFlVvZrk88BB4H+Bx4C9wBtVdaTbbBZYvdjzk2wBtgBccsklo46hU8wPF5OWnz6Xbs4HNgCXAu8CzgFuWGTTWuz5VbW9qmaqamZqamrUMSRJx9Hn0s0Hge9V1VxV/Qz4CvB+4LzuUg7AGuC1njNKknroE/qDwDVJzk4SYB3wIvAE8JFum03Arn4jSpL6GDn0VfUk8990fQp4rnut7cBtwGeSHADeCdw7hjklSSPq9RumquoO4I6jVr8EXN3ndSVJ4+M7YyWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcb1Cn+S8JDuTfCfJviTvS3JBkseT7O9uzx/XsJKkk9f3jP6LwNeq6j3Ae4F9wFZgd1WtBXZ3y5KkgYwc+iTnAh8A7gWoqjer6g1gA7Cj22wHsLHvkJKk0fU5o383MAd8KcnTSe5Jcg5wUVUdAuhuLxzDnJKkEfUJ/QrgKuDuqroS+AkncZkmyZYke5LsmZub6zGGJOnt9An9LDBbVU92yzuZD//rSVYBdLeHF3tyVW2vqpmqmpmamuoxhiTp7Ywc+qr6PvBKksu6VeuAF4GHgU3duk3Arl4TSpJ6WdHz+Z8EHkhyJvAS8HHmv3g8lGQzcBC4uec+JEk99Ap9VT0DzCzy0Lo+rytJGh/fGStJjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9Jjev7oWYawPTWR4YeQdIy4hm9JDXOM3pJS8ZQf1t9edv6QfZ7qnhGL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1LjeoU9yRpKnk3y1W740yZNJ9if5cpIz+48pSRrVOM7obwX2LVi+E7irqtYCPwI2j2EfkqQR9Qp9kjXAeuCebjnAdcDObpMdwMY++5Ak9dP3jP4LwGeBn3fL7wTeqKoj3fIssHqxJybZkmRPkj1zc3M9x5AkHcvIoU9yI3C4qvYuXL3IprXY86tqe1XNVNXM1NTUqGNIko6jz++MvRa4KcmHgbOAc5k/wz8vyYrurH4N8Fr/MSVJoxr5jL6qbq+qNVU1DdwCfL2qPgY8AXyk22wTsKv3lJKkkU3i5+hvAz6T5ADz1+zvncA+JEknqM+lm/9XVd8AvtHdfwm4ehyvK0nqz3fGSlLjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNW7F0AMsZ9NbHxl6BEk6Ls/oJalxhl6SGmfoJalxhl6SGjdy6JNcnOSJJPuSvJDk1m79BUkeT7K/uz1/fONKkk5WnzP6I8AfV9VvAtcAn0hyObAV2F1Va4Hd3bIkaSAjh76qDlXVU939/wb2AauBDcCObrMdwMa+Q0qSRjeWa/RJpoErgSeBi6rqEMx/MQAuHMc+JEmj6R36JL8C/CPwR1X1XyfxvC1J9iTZMzc313cMSdIx9Ap9kl9mPvIPVNVXutWvJ1nVPb4KOLzYc6tqe1XNVNXM1NRUnzEkSW+jz0/dBLgX2FdVf7ngoYeBTd39TcCu0ceTJPXV57NurgV+D3guyTPduj8FtgEPJdkMHARu7jeiJKmPkUNfVf8K5BgPrxv1dSVJ4+U7YyWpcYZekhrn59FLOu0N+bslXt62fuL78Ixekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpccv+IxCGfOuyJC0HntFLUuMMvSQ1ztBLUuMMvSQ1ztBLUuMMvSQ1ztBLUuMMvSQ1ztBLUuMMvSQ1ztBLUuMMvSQ1ztBLUuMMvSQ1biKhT3J9ku8mOZBk6yT2IUk6MWMPfZIzgL8BbgAuBz6a5PJx70eSdGImcUZ/NXCgql6qqjeBfwA2TGA/kqQTMInQrwZeWbA8262TJA1gEr9KMIusq7dslGwBtnSL/5Pku0dtshL4wZhnW6pOl2M9XY4TPNYWTeQ4c2evp//6iWw0idDPAhcvWF4DvHb0RlW1Hdh+rBdJsqeqZsY/3tJzuhzr6XKc4LG2aDkf5yQu3fw7sDbJpUnOBG4BHp7AfiRJJ2DsZ/RVdSTJHwL/ApwB3FdVL4x7P5KkEzOJSzdU1aPAoz1f5piXdRp0uhzr6XKc4LG2aNkeZ6re8n1SSVJD/AgESWrckgx9kjOSPJ3kq0PPMklJzkuyM8l3kuxL8r6hZ5qUJJ9O8kKS55M8mOSsoWcalyT3JTmc5PkF6y5I8niS/d3t+UPOOA7HOM6/6P78Ppvkn5KcN+SM47LYsS547E+SVJKVQ8w2iiUZeuBWYN/QQ5wCXwS+VlXvAd5Lo8ecZDXwKWCmqq5g/pv0tww71VjdD1x/1LqtwO6qWgvs7paXu/t563E+DlxRVb8F/Adw+6keakLu563HSpKLgQ8BB0/1QH0sudAnWQOsB+4ZepZJSnIu8AHgXoCqerOq3hh2qolaAbwjyQrgbBZ5b8VyVVXfBH541OoNwI7u/g5g4ykdagIWO86qeqyqjnSL32L+fTPL3jH+mwLcBXyWRd4EupQtudADX2D+X+TPhx5kwt4NzAFf6i5T3ZPknKGHmoSqehX4PPNnQYeAH1fVY8NONXEXVdUhgO72woHnORX+APjnoYeYlCQ3Aa9W1beHnuVkLanQJ7kROFxVe4ee5RRYAVwF3F1VVwI/oY2/3r9Fd316A3Ap8C7gnCS/O+xUGqcknwOOAA8MPcskJDkb+BzwZ0PPMoolFXrgWuCmJC8z/6mX1yX5+2FHmphZYLaqnuyWdzIf/hZ9EPheVc1V1c+ArwDvH3imSXs9ySqA7vbwwPNMTJJNwI3Ax6rdn9f+DeZPVL7d9WkN8FSSXxt0qhO0pEJfVbdX1Zqqmmb+m3Vfr6omz/yq6vvAK0ku61atA14ccKRJOghck+TsJGH+WJv8xvMCDwObuvubgF0DzjIxSa4HbgNuqqqfDj3PpFTVc1V1YVVNd32aBa7q/j9e8pZU6E9DnwQeSPIs8NvAnw88z0R0f2vZCTwFPMf8n7tl+y7DoyV5EPg34LIks0k2A9uADyXZz/xPaWwbcsZxOMZx/jXwq8DjSZ5J8reDDjkmxzjWZct3xkpS4zyjl6TGGXpJapyhl6TGGXpJapyhl6TGGXpJapyhl6TGGXpJatz/AbDD8Ko/TeGRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(np.matmul(A,np.ones((n,1))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
